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■ ABSTRACT 

ir^ ' Rigorous computer simulations of propagating electromagnetic fields have become an important tool for optical 
CIh metrology and optics design of nanostructured components. As has been shown in previous benchmarks some of 
the presently used methods suffer from low convergence rates and/or low accuracy of the results and exhibit very 
: long computation times^' ^ which makes application to extended 2D layout patterns impractical. We address 3D 
^ ■ simulation tasks by using a finite-element solver which has been shown to be superior to competing methods by 
, several orders of magnitude in accuracy and computational time for typical microlithography simulations.^ We 
report on the current status of the solver, incorporating higher order edge elements, adaptive refinement methods, 
; and fast solution algorithms. Further, we investigate the performance of the solver in the 3D simulation project 
• ' of light diffraction off an alternating phase-shift contact-hole mask. 

in ; 
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1. INTRODUCTION 

With the advances of micro- and nanotechnology simulation tools for rigorous solutions of Maxwell's equations 
have become an important tool in research and development. Designing, e.g., a nanooptical component or a 
d metrology tool is usually assisted by computer simulations, and simulations are used in most scientific works on 
nanooptical research to support theoretical and experimental findings. It exists a variety of different methods for 
solving Maxwell's equations, and generally also a variety of different numerical implementations of each method. 
Prominent examples of different methods are the finite-element method (FEM) , the finite difference time domain 
method (FDTD), the boundary element method (BEM), and rigorously coupled wave analysis (RCWA). 

We address 3D simulation tasks occuring in microlithography by using the frequency-domain finite-clement 
solver JCMsuite. This solver has been successfully applied to a wide range of 3D electromagnetic field compu- 
tations including microlithography,^"^ left-handed metamaterials in the optical regime,^' ^ photonic crystals,^ 
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and ncarficld-microsc:opy7 The solver has also been used for pattern reeonstruetion in EUV scatterometry,^' ® 
and it has been benehniarked to other methods in typical DUV lithography^'^ and other*"' projects. 

In this paper we report on the current status of the finite-element solver JCMsuite, and we present simulations 
of light transition through periodic arrays of alternating phase-shift contact-holes which have to be printed on 
wafers with the same reticle with high overlapping process window. 3D effects cause imbalances of intensities of 
the 0°- and 180°-holes, and of different behavior in defocus if the etch depth is not adjusted. The evaluation of 
the impact of these effects has been chosen as a test case for JCMsuite. 

2. BACKGROUND 

The finite-element package JCMsuite allows to simulate a variety of electromagnetic problems. It incorporates a 
scattering solver (JCMharmony) , a propagating mode solver (JCMmode) and a resonance solver (JCMresonance) . 
The scattering, eigenmode and resonance problems can be formulated on ID. 2D and 3D computational domains. 
Admissible geometries can consist of periodic or isolated patterns, or a mixture of both. Further, solvers for 
problems posed on cylindrically symmetric geometries are implemented. 
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Figure 1. Schematics of the 3D test structure: cross-section in a x-y-plane. 0-degree (180-degree) rectangular phase-shift 
holes are depicted in grey (light grey). The elementary vectors, ai,02, of the periodic pattern are indicated by arrows. 

In this paper we concentrate on light-scattering off a 3D pattern (photomask) which is periodic in the x— 
and y— directions and is enclosed by homogeneous substrate (at Zsut) and superstrate (at Zsup) which arc infinite 
in the — , resp. H-z— direction (see Figures 1 and 2). Light propagation in the investigated system is governed 
by Maxwell's equations where vanishing densities of free charges and currents are assumed. The dielectric 
coefficient e(x) and the permeability /i(x) of the considered photomasks are periodic and complex, s (x) = 
e {x + a), jjL (x) = ii{x + a). Here a is any elementary vector of the periodic lattice. For given primitive lattice 
vectors Si and a2 an elementary cell C can be defined as O = {x G | a; = aiUi + 0202; < ai, 0:2 < l} x 
[zsub,Zsup]- As can be seen from the computational domain in Figure 3 elementary cells 12 of the same volume 
can be defined also by more complicated polygons which are adapted to the actual geometry. * 

*In this case a computational domain constructed corresponding to Q, = {x ^ \ x = aiSi + 02(12; < Qi, Q2 < l} x 
[zsub,Zeup] would in general lead to intersections of the domain boundary and geometrical features at sharp angles. This 
could result in low quality triangulations. 
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Figure 2. Schematics of the 3D test structure: cross-section in a x-z-planc. The material stack on a Quartz (Si02) 
substrate consists of a molydenum sihcide and two layers containing chromium. 0-degree (180-degree) rectangular phase- 
shift air-filled holes are indicated. 
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Figure 3. Schematics of the 3D test structure: cross-section in a a;-«-plane. The computational domain is indicated by 
a polygon. 

A timc-harmonic ansatz with frequency w and magnetic field H{x,t) = e~^^H{x) leads to the following 
equations for H{x): 



The wave equation and the divergence condition for the magnetic field: 

V X ^— V X H(x) - J^n{x)H{x) = 0, f e O, 
e{x) 

V • fi{x)H{x) = 0, xGfl. 



(1) 
(2) 



Transparent boundary conditions at the boundaries to the substrate (at Zgub) and superstrate (at z^up), 
dCl, where iJ'" is the incident magnetic field (plane wave in this case), and n is the normal vector on dfl: 



s{x) 



V X {H - W'')] xn = DtN{H - W''), x G dft. 



(3) 



The DtN operator (Dirichlet-to-Neumann) is realized with an adaptive PML method. ^'^ This is a 
generalized formulation of Sommerfeld's radiation condition. 



• Periodic boundary conditions for the transverse boundaries, dfl, governed by Bloch's theorem^^: 

H{x) = e'^-^w(f), u{x) = u{x + a), (4) 

where the Bloch wavcvcctor fc G R"^ is defined by the incoming plane wave if™. 

Similar equations arc foimd for the electric field E(x.t) = e^'^'^^ E{x)\ these are treated accordingly. The 
finite-element method solves Eqs. (1) - (4) in their weak form, i.e., in an integral representation. 

The finite-element methods consists of the following steps: 

• The computational domain is discretized with simple geometrical patches, JCMsuite uses linear (ID), 
triangular (2D) and tetrahedral or prismatoidal (3D) patches. The use of prismatoidal patches is well 
suited for layered geometries, as in photomask simulations. 

• The function spaces in the integral representation of Maxwell's equations are discretized using Ned- 
elec's edge elements, which are vectorial functions of polynomial order defined on the simple geometrical 
patches. In the current implementation, JCMsuite uses polynomials of first (!'' ) to ninth (9*^*) order. In 
a nutshell, FEM can be explained as expanding the field corresponding to the exact solution of Equation (1) 
in the basis given by these elements. 

• This expansion leads to a large sparse matrix equation (algebraic problem) . To solve the algebraic problem 
on a standard workstation linear algebra decomposition techniques (LU- factorization, e.g.. package PAR- 
DISO,^^ which was used in the simulations of Chapter 3) are used. In cases with either large computational 
domains or high accuracy demands, also domain decomposition methods^'' are used and allow to handle 
problems with very large numbers of unknowns. 

For details on the weak formulation, the choice of Bloch-periodic functional spaces, the FEM discretization, 
and the implementation of the adaptive PML method in JCMsuite we refer to previous works. In future im- 
plementations performance will further be increased by using curvilinear elements, general domain-decomposition 
techniques and /ip- adaptive methods. 

3. 3D SIMULATIONS OF A CONTACT-HOLE PHOTOMASK 
3.1. Parameter settings and simulation flow 

We apply the 3D light-scattering solver of the programme package JCMsuite to the test case of a 2D periodic 
pattern of an alternating phase-shift contact-hole photomask. Figures 1, 2 and 3 show cross-sections through the 
geometrical layout of the mask and define the geometrical parameters. The geometrical parameters investigated 
in this study are given in Table 1. The computational domain is a unit-cell of the periodic pattern; it is indicated 
in Figure 3. Its total volume is = 2px x py x h^, where = cIq^. + cIq-^ bottom "^MoSi ^~ "^etch" -Please 
note that a simple rectangular box- like computational domain would have twice the volume, V = 2px x 2py x hz- 
Within the numerical error the results are independent of the choice of the unit-cell chosen as computational 
domain. 

Figure 4 shows the simulation flow of the finite-element software: 

• The geometry of the computational domain is described in a polygonal format of the cross-section, cor- 
responding to Fig. 1, and includes a height profile (see Fig. 2) with material attributions and statements 
about the computational domain boundaries (periodic and transparent boundaries in this case). The 
translational vectors of the periodic pattern (ai,a2) are identified automatically from the layout, opti- 
mized settings of the perfectly matched layers (PML) are found automatically with an adaptive method 
(adaptive PML, aPML). Further, parameters specifying a maximum patch size of the finite-element mesh 
and further meshing properties can be set. From these geometry parameters and mesh parameters the pris- 
matoidal mesh is generated automatically. Figure 5 shows {x — t/)-cross-sections through the prismatoidal 
meshes for different geometrical parameter sets. 
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Figure 4. Schematics of tlie simulation flow of the FEM solver JCMsuite. 
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set 1 


set 2 


set 3 


set 4 


set 5 


set 6 


Px [nm] 
Py [nm] 
CDx,0 [nm] 
CDy [nm] 


576 
1120 
416 
876 


640 
1120 
452 
800 


720 
1120 
496 
752 


800 
1120 
552 
696 


1040 
1120 
496 
976 


1200 
1120 
520 
980 


CDx,l80 [nm] 
CDy, 180 [nm] 
3D bias [nm] 

dCr,top [nm] 

dCr,bottom [nm] 
^MoSi [nm] 
^ctch [mn] 


CDx,0 + 3D bias 

CDy + 3D bias 

(32^ 48, 64, 80) 

20 

40 

68 

(V,(^, 160, 164, 168, 172, 176, 180] 


Ao 

^rair 
^rCr,top 
^rCr,bottom 
SrMoSi 

£rSi02 


193.0 nm 
1.0 

(1.871 + 1.13i)2 
(1.477+ 1.762i)2 
(2.343 + 0.586i)2 
(1.564)2 



Table 1. Parameter settings for the 3D simulations. The relative permittivity, er, is the square of the complex index of 
refraction, Sr = = {rir +ik)'^. The simulation wavelength was A = 193.0 nm while the wavelength in the real system is 
A = (193.36 ± 0.2) nm. The dispersion of e is taken into account at A = 193.36 nm. All combinations of 3D bias and etch 
depth, dg^j,jj, have been investigated for the six parameter sets in order to find for each pattern the best value and the 
best overall 3D bias and etch depth. 



• Material parameters (complex permittivity and permeability tensors) can be specified as piecewise constant 

functions and/or (using dynamically loaded libraries) as functions of arbitrary spatial dependence. In the 
presented example, the piecewise constant, isotropic settings given in Table 1 are used. The relative 
permeability is /i^ = 1 for all used materials. 

• Light sources can be defined as predefined functions (plane waves, Gaussian beams, point sources) or as 

arbitrary functions using dynamically loaded libraries. JCMsiiite allows to generate solutions to several 
independent source terms in parallel, efficiently re-using the inverted system matrix. In the presented exam- 
ple, two plane waves with normal incidence (wavevector k = (0, 0, fc^)), and with orthogonal polarizations 
are used. 

• The main project definitions in this case are the accuracy settings (mesh refinement, PML refinement, 
finite-element degree) and the definitions of postprocesses. Upon execution the JCMsolve computes the 
full field distribution over the entire computational domain. Through internal or external post-processes, 




Figure 5. Cross-sections through prismatoidal discretizations of several computational domains with different geometry 
parameters, corresponding to parameter sets 1 (a), 4 (b), 6 (c)(see Table 1). Please note that (a), (b), and (c) have 
different length scales. 

the quantities of interest can be derived from the field. E.g., the complex amplitudes of propagating modes 
are attained by Fourier transformation, an aerial image is calculated,'^ or the field distribution is exported 
to graphics format for visualisation and analysis. 

• Interfaces to scripting languages like Python or MATLAB can be used for performing automatic parameter 
scans and data analysis. In the presented example, the MATLAB interface for the generation of the input 
files and for the execution of a loop over the geometrical parameters given in Table 1 consists of well below 
100 lines of MATLAB code. 

3.2. Error analysis 

We have investigated the accuracy of the near field results by simulating near fields using finite-elements of differ- 
ent polynomial degrees (p-refinement) and using meshes with different numbers of mesh elements (/i-refinement) 
for one fixed setting of geometrical and other physical parameters. From each near fields we compute the diffrac- 
tion pattern using a post-process. Figure 6 shows the convergence of the absolute values of the intensities in the 
zero and first diffraction orders. Figure 7(a) shows the convergence of the relative error of the intensity / in the 
5x5 lowest diffraction orders. The relative error of the intensity, A/^, is defined as (/ — Iqe)/Iqe, where the 
quasi-exact intensity, Iqe , is deduced from a solution with highest finite-element degree and finest mesh. As can 
be seen from Figure 7, relative errors of below one percent for the prominent diffraction orders can be reached 
using a high polynomial degree of the finite-elements {p > 4) and relatively moderate numbers of unknowns, 
corresponding to coarse spatial discretizations. 

3.3. Computational costs 

Figure 7(b) shows the computation time in minutes in dependence on the number of degrees of freedom in 
the simulation. Different symbols in the plot correspond to different finite-element polynomial degrees. The 
computations have been performed on a standard workstation with 16 (2 x 8) processors and extended RAM 
(AMD Opteron, 64GB RAM). Please note that for the results in Figure 7(b) we have used only four of the 
16 processors, while for the scans over geometrical parameters we have used twelve out of the 16 processors, 
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Figure 6. Convergence of the intensities of the first (a) and zero (b) diffraction orders with increasing number of unknowns 
of the FEM solution. Results from computations with finite element degree p = 2 . . .6 are displayed. 
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Figure 7. (a) Convergence of the total transmission amplitude to the 5x5 lowest diffraction orders, (b) Computation 
times in dependence of the number of unknowns of the corresponding matrix equation. 



giving a speed-up factor of approximately three in the computation times. As can be seen in Figure 6 and 
Figure 7, results corresponding to an accuracy of the prominent diffraction orders of better 1% can be attained 
in computation times of few minutes when using finite elements of order p = 4 and above. Memory requirements 
were roughly 5. . . 15GB RAM for A'' ~ 2.5 ... 5 • 10^ and finite-element degree p = 4 (the setting used for the scan 
over the geometrical parameters) . Please note that in general it does depend on the physical problem wheather 
a refinement of the finite element degree p, p p + I, or & refinement of the mesh width h, h ^ h/2, leads 
to a better convergence of the results. For the given rather large (on the wavelength-scale) structures on the 
photomask, rather high degrees p seem to lead to the best results, while for problems with smaller features finer 
meshes are favourable. In a future implementation of the programme package this choice will be automatically 
and locally detected by /ip- adaptive methods. 
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Figure 8. Intensity of the first (+1,0) diffraction order for the different investigated parameter sets. Each 
parameter set consists of four groups (3D bias=8, 12, 16,20 nm) of seven values of etch depths (dg^^jj = 
156, 160, 164, 168, 172, 176, 180 nm). 
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Figure 9. Intensity of the zero (0,0) diffraction order for the different investigated parameter sets. Each parameter set con- 
sists of four groups (3D bias=8, 12, 16,20 nm) of seven values of etch depths (dg^^jj = 156, 160, 164, 168, 172, 176, 180 nm). 

3.4. Geometry dependence of the diffraction spectrum 

We have performed a scan over the geometrical parameters defined in Table 1. The six parameter sets consist 
each of 28 subsets of the various combinations of 3D bias and etch depth. Figure 8 and 9 show how the first 
(+1,0) order and zero order diffraction intensities vary over the scan of geometrical parameters. As expected for 
phase-shift masks, the intensity in the first diffraction orders (prominent order) is about two orders of magnitude 
higher than the intensity in the zero diffraction order. For otherwise fixed geometrical parameters, the variation 
of the prominent diffraction order with the etch depth in the investigated regime is about 2-4%. The variation 
of the prominent diffraction order with the 3D bias is about 10-15%. The variation with a different set of CD 
and pitch parameters is much larger. 

This shows that for mask design using rigorous simulations, especially for the design of parameters such as 
the etch depth, it is necessary to use simulation tools which produce results of low numerical error. 
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Figure 10. (a) CD as a function of defocus for an almost balanced case, (b) CD as a function of defocus for an imbalanced 
case. 



3.5. Imaging application 

For the conditions of Table 1 the farfield coefRcients were determined and used as input for an Qimonda inhouse 
imaging simulator. The mask bias for each pitch was pre-optimized by Kirchhoff (2D) mask simulations in 
resist. For evaluation of the arial image output for each pitch and 3D bias the intensity threshold for target CD 
was determined at best focus as average of 0°- and 180°-holes. Figure 10(b) shows the typical phase balancing 
problem that occurs in case of non-optimum etch depth in the mask: the CD of 0°- and 180°-holes behaves 
different at defocus, whereas for the optimum depth the behavior is symmetric (Figure 10(a)). 

The CD difference of 0°- and 180°-holes over defocus shows a slope that is characteristic for the phase error, 
i.e., the deviation from optimum etch depth where the defocus behavior is flat (see Figure 11(a)). Furtheron, the 
180°-holc prints smaller than a 0°-holc of same size due to 3D mask effects. This is compensated by a 3D mask 
bias, see Figures 11(b) and 12(a), which is in the range of 8-20 nm, depending on pattern and mask conditions 
(side wall angles, material). The slope is rather independent of the 3D bias, as shown in Figure 12(b). The 
optimum 3D bias, where the CD difference is zero, is 8.5 nm at an etch depth of 168 nm, as obtained by a linear 
fit of the CD difference vs. 3D bias. As can be also seen in Figure 12(b), the slope is almost zero at this etch 
depth. 

A determination of the optimum etch depths and optimum 3D mask biases through pitch reveals significant 

differences. While the 3D mask bias can be corrected individually for such special patterns by software algorithms 
it is challenging to correct the 3D mask effect in real layouts containing random patterns. An even more critical 
issue is to find an optimum etch depth through pitch. Here an analysis of overlapping process windows including 
etch depth deviations is the next step to find an optimum etch depth through pitch. Then conclusions can be 
drawn with respect to lithographic applications for nodes below 40 nm half pitch. 



4. CONCLUSIONS 

We have performed rigorous 3D FEM simulations of light transition through alternating phase-shift contact- 
hole photomasks using the finite-element programme package JCMsuite. We have investigated the convergence 
behavior of the solutions and we have shown that we achieve results at high numerical accuracy. Our results 
show that rigorous 3D mask simulations can well be handled at high accuracy and relatively low computational 
cost. 

The investigation of the diffraction spectrum of an alternating phase-shift contact-hole mask has shown that 
for mask design and optimization projects it is necessary to use rigorous tools with low discretization errors, 
typically below 1% for the prominent diffraction orders, because the variance of the diffraction orders can be 
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Figure 11. (a) CD difference (CDq — CDj^gg) as a function of defocus and etch depth for one fixed 3D bias, (b) CD 
difference (CDq — CDj^gg) as a function of defocus and 3D bias at optimized etch depth. 
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Figure 12. (a) CD difference (CDq — CDj^gg) as a function of 3D bias at optimized etch depth (168 nm). (b) Slope of 
CD difference, 9(CDq — CD]^gQ)/9(defocus), as a function of etch depth and 3D bias. 



of the order of few percent for typical design parameter regimes. We have used the simulations to determine 
optimum etch depths and 3D mask biases for the various investigated parameter sets. 
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